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INTRODUCTION 


Modal  techniques  are  frequently  employed  to  solve  the  free  vibrational 
motion  of  beams.  For  uniform  beams,  the  method  is  analytical  and  is  described 
in  the  textbook  literature  (refs  1-4).  However,  when  the  beam's  cross  section 
is  nonuniform,  it  is  often  difficult,  if  not  impossible,  to  find  an  analytical 
solution.  One  of  the  most  commonly  used  methods  for  beams  in  this  category  is 
to  divide  the  beam  into  a  finite  number  of  uniform  sections  and  solve  the 
equations  within  each  section.  Klein  (ref  5)  developed  a  "component  inodes" 
method  based  on  this  concept  which  she  used  to  analyze  a  cantilevered  helicopter 
blade  with  a  high  degree  of  nonuniformity.  For  each  section  of  the  blade,  modal 
displacements  were  expressed  by  a  Rayleigh-Ritz  expansion,  and  continuity 
between  sections  was  enforced  by  Lagrange  multipliers.  Mittendorf  and  Greif 
(ref  6)  used  a  similar  approach  for  a  beam  with  general  boundary  conditions. 

The  mode  shapes  were  expanded  in  a  Fourier  series,  while  the  Lagrange 
multipliers  were  used  to  enforce  the  kinematic  transfer  across  section  bound¬ 
aries.  A  finite  element  approach,  such  as  the  one  developed  by  Resende  and 
Doyle  (ref  7),  may  also  be  used  for  these  cases.  They  used  a  three-noded  ele¬ 
ment  and  a  formulation  based  upon  the  quadratic  isoparametric  class.  Identical 
shape  functions  are  used  to  define  geometry,  displacement,  and  variations  in 
cross-secti-onal  properties. 

For  beams  with  special  geometrical  properties,  a  number  of  analytical  tech¬ 
niques  may  be  used.  Wang  (ref  8)  obtained  a  general  solution  for  beams  whose 
sectional  properties  vary  as  a  power  of  the  axial  coordinate.  Mode  shapes  were 
obtained  in  terms  of  generalized  hypergeometric  functions  using  the  method  of 
Frobenius.  Wang  (ref  9)  used  the  Galerkin  method  to  determine  frequencies. 

Modes  were  represented  by  a  set  of  Legendre  polynomials  which  served  as  shape 
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functions  throughout  the  spatial  domain.  Lau  (ref  10)  derived  the  frequency 
determinant  for  a  profiled  and  truncated  fixed-free  beam  with  an  attached  end 
mass  using  Bessel  functions  of  the  second  order.  Goel  (ref  11)  investigated 
vibrations  of  linearly-tapered  beams  mounted  upon  rotationally-elastic  bound¬ 
aries. 

All  of  these  methods  are  satisfactory  within  the  confines  of  their 
modelling  approximations.  For  example,  with  finite  elements,  the  mode  shape 
between  nodes  is  approximated  by  a  piecewise  set  of  cubic  polynomials. 

Continuity  at  element  boundaries  is  satisfied  up  to  the  second  derivative,  and 
in  most  applications  this  is  all  that  is  required.  Unlike  the  FEM,  analytical 
techniques  use  functional  representations  which  are  particularly  useful  in 
applications  where  differentiability  is  required. 

The  purpose  of  this  work  is  to  present  a  hybrid  method  (analytical/ 
numerical)  to  model  the  free  vibrations  of  statically-determinant  beams  on 
general  supports.  This  method  adapts  the  exact  prismatic  beam  solution  to 
nonprismatic  cases.  In  addition  to  a  presentation  of  this  method,  the  computer 
routines  that  implement  the  method  are  explained  and  a  few  examples  are  compared 
to  the  results  generated  by  the  FEM  method. 


ANALYTICAL  METHOD 

The  transverse  motion  of  a  general  nonprismatic  beam  is  described  by  the 
following  partial  differential  equation  (PDE)  according  to  Thomson  (ref  1): 
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where 


E  =  Young's  Modulus 

I  =  transverse  moment  of  inertia 

J  =  pitch  moment  of  inertia 

G  «  modulus  of  elasticity  in  shear 

A  =  area  of  beam's  cross  section 

k  =  shear  coefficient  (function  of  cross  section) 

w  =  weight  per  unit  length  of  beam 

pn  a  n-th  forcing  function  (total  of  N) 

g  =  acceleration  of  gravity 

x  =  axial  coordinate 

t  -  time  (independent  variable) 

y  =  transverse  displacement  (dependent  variable) 

•  =  time  derivative 

'  =  space  derivative 

This  equation  takes  into  account  all  possible  effects  of  beam  dynamics.  The 
second  term  accounts  for  the  translational  inertia  of  the  beam,  while  the 
remaining  terms  model  rotatory  inertia  and  shear  deformation  effects.  All  terms 
are  on  a  per  length  basis.  The  right  side  contains  the  driving  loads  which  may 
be  functions  of  the  dependent  variable  as  well  as  space  and  time. 

The  free  vibration  modes  (or  mode  shapes)  are  determined  by  solving  the 
homogeneous  form  of  Eq.  (1).  For  a  prismatic  structure  (i.e.,  when  I,  E,  J,  G, 
A,  w,  k  are  constant  throughout  the  spatial  domain),  the  free  vibration  equation 
takes  the  form 
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The  solution  of  this  equation  is  straightforward  and  can  be  found  in  the  works 
of  Bozich  (ref  12)  and  Dolph  (ref  13).  Following  the  method  used  by  Bozich,  the 
POE  becomes  a  function  of  the  structural  parameters  and  the  natural  vibration 
frequencies  as  follows: 
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where 

w  =  natural  vibration  frequency 

Equation  (3)  can  be  nondimensional ized  by  introducing  a  new  independent  variable 
x  (=x/L).  By  normalizing  each  term  with  respect  to  w2/g  and  using  the  parame¬ 
ters  suggested  by  Kruszewski  (ref  14),  the  resulting  dimensionless  equation, 
valid  for  any  uniform  beam,  is 
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These  coefficients  represent  the  contributions  of  bending,  shear,  and  rotatory 
inertia,  respectively.  In  this  form,  Eq.  (4)  has  the  solution: 

y(x,t)  =  [A  cosh  j3x  +  B  sinh  0x  +  C  cos  ax  +  D  sin  ax]cos  ut  (6) 


where 
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The  coefficients  A,  B,  C,  and  D  are  determined  by  the  boundary  conditions  of  the 
beam.  The  vibration  frequencies,  w,  and  their  respective  mode  shape  functions 
are  found  by  applying  the  appropriate  boundary  conditions. 

To  adapt  this  analysis  to  the  nonprismatic  case,  the  beam  is  represented  by 
a  finite  number  of  prismatic  segments.  Within  each  segment  Eq.  (6)  portrays  the 
mode  shape  for  the  assumed  uniform  properties  of  that  segment.  The  conditions 
at  the  common  boundary  between  adjacent  segments  must  match  to  preserve  con¬ 
tinuity  between  segments.  The  appropriate  end  boundary  conditions  are  applied 
at  the  extremities  of  the  beam.  The  equation  defining  the  mode  shapes  for  each 
segment  is  developed  from  Eq.  (6)  as  follows-. 

d> T  j  ( x )  a  A-j  j  cosh  0-jjX  +  B-jj  sinh  /3-j  jx  +  C-jj  cos  a-jjx  +  D-jj  sin  a-jjx  (8) 

where 

i  =  segment  number 
j  =  mode  shape  number 

A , B , C , D  =  function  coefficients  (by  mode/segment) 

Employing  these  segmented-mode  shape  functions,  the  natural  vibration  frequen¬ 
cies  and  mode  shape  coefficients  can  be  determined  to  within  an  arbitrary  con¬ 
stant  for  any  set  of  end  boundary  conditions.  For  example,  consider  the  case  of 
a  free-free  beam  which  is  composed  of  M  segments.  The  boundary  conditions  are 
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Continuity  at  the  common  boundaries  between  the  segments  is  preserved  by 
equating  displacement,  slope,  moment,  and  shear  calculated  by  adjacent  mode 
functions.  These  conditions  are  as  follows: 

4>i-l j ( x i )  =  4>ij(Xi)  (10a) 

4>i  —  ij(xi)  =  i  j  ( x  t  )  (10b) 

(El )  i  —  1  *l>i  —  1  j  ( x i )  =  (El )  i4>i  j  (x-j )  (10c) 

(EDi.^-lTjfXi)  =  (EDi^jfxi)  (lOd) 

where 

x^  =  left  boundary  of  i-th  segment 

For  a  model  containing  M  segments,  there  are  4M  coefficients  and  one  frequency 
to  be  evaluated  at  each  mode  shape.  Setting  a  coefficient  of  the  mode  shape  in 
the  first  segment  to  unity  renders  the  system  deterministic  and  yields  unique 
solutions  at  each  natural  frequency. 

A  set  of  algebraic  equations  in  terms  of  the  natural  frequency,  ,  and  the 
4M-1  coefficients  is  presented  in  matrix  form  in  Figure  1.  The  system  matrix 
contains  an  orderly  set  of  entries  representing  the  enforcement  of  the  boundary 
conditions  (B.C.)  and  continuity  conditions  (C.C.).  At  global  location  1,1  the 
2  by  4  subarray  contains  the  terms  for  evaluating  the  B.C.  at  x  =  0,  whereas  a 
similar  array  ending  at  N,N  represents  the  B.C.  at  x  =  1.  "Walking"  along  the 
main  diagonal  are  M-l  subarrays,  4  by  8  in  size,  which  represent  the  C.C.  across 
segment  boundaries.  The  global  address  of  the  upper  left  corner  of  the  initial 
subarray  is  3,1  and  is  incremented  4,4  for  each  additional  interface.  The 
remaining  terms  are  zero.  The  nonzero  entries  are  functions  of  the  segment 
properties,  end-point  locations,  and  the  vibration  frequencies.  Setting  the 
determinant  equal  to  zero  and  solving  for  the  roots  of  the  resulting  charac¬ 
teristic  equation  produces  values  for  these  frequencies.  At  each  frequency,  a 
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reduced  system  is  developed  by  setting  the  appropriate  mode  shape  coefficient  in 
the  first  segment  to  unity,  eliminating  the  first  row,  and  shifting  the 
corresponding  column  to  the  right  side  of  the  matrix  equation.  This  new  system 
of  equations  is  deterministic  at  each  frequency.  The  values  for  the  remaining 
mode  shape  coefficients  may  then  be  calculated.  The  orderliness  of  the  entries 
to  the  system  matrix  enables  an  efficient  and  dedicated  computer  code  for  its 
solution  to  be  developed. 

This  method  has  one  main  advantage  over  finite  element  techniques  because 
the  segmenting  of  naturally  uniform  sections  is  not  required  since  the  mode 
function  is  exact  within  this  domain.  Thus,  only  a  few  segments  are  needed  for 
structures  with  uniform  or  gradually  varying  sections.  Another  advantage  of 
this  method  is  that  the  number  of  modes  which  can  be  calculated  is  not  a  func¬ 
tion  of  the  numher  of  segments.  In  finite  elements,  the  mode  shapes  are  defined 
by  the  location  of  the  model's  nodes  so  that  an  accurate  representation  of  the 
mode  shapes  at  higher  frequencies  requires  a  large  number  of  nodes.  This  is  not 
the  case  for  the  USM,  since  its  "nodeless"  elements  use  trigonometric  and  hyper¬ 
bolic  functions  to  represent  the  mode  shapes.  Finally,  in  cases  where  the  modes 
are  used  in  dynamic  analyses  (ref  15)  and  the  loads  are  functions  of  the  depend¬ 
ent  variable  and  its  derivatives,  the  USM  mode  functions  can  accurately  repre¬ 
sent  these  loads  since  they  are  differentiable.  The  orthogonality  of  the  USM- 
generated  mode  shapes  permits  the  use  of  the  usual  expansion  techniques.  Their 
finite  element  counterparts  are  usually  represented  by  cubic  polynomials  which 
possess  limited  differentiability. 

COMPUTER  CODES 

A  set  of  FORTRAN-77  computer  routines  has  been  written  and  tested  on  an  IBM 
4341-N12  minicomputer  driven  by  the  VM/SP  Operating  System.  A  total  of  2000 
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lines  of  code  was  written  for  this  analysis.  The  overall  operation  is 
controlled  by  an  executive  routine  called  MOD-MAIN.  This  routine  directly  calls 
7  of  the  20  subroutines  used  in  the  analysis.  Figure  2  is  a  flowchart  of  the 
overall  operation.  The  data  entry  and  initialization  portion  of  the  analysis  is 
shown  on  the  left  side  of  the  figure.  Specifically,  input  data  consists  of  the 
physical  and  geometric  characteristics  of  the  beam,  the  location  of  segment  bor¬ 
ders,  and  the  boundary  conditions  at  the  beam  extremities.  In  addition,  the 
number  of  modes  to  be  determined,  parameters  regarding  the  search  procedure, 
tolerance  on  the  root  estimates,  and  output  requirements  are  read  from  a  stored 
disk  file. 

The  right  section  of  Figure  2  shows  the  flow  through  the  analytical  portion 
of  the  routines.  A  loop  is  shown  that  is  traversed  once  for  each  mode 
requested.  During  each  pass,  a  vibration  frequency  and  set  of  mode  shape  coef¬ 
ficients  are  determined  and  the  results  sent  to  various  disk  files.  The  proce¬ 
dure  for  locating  a  root  is  contained  within  the  routine  labelled  ROOT  SOLVER. 

In  this  solver,  a  search  procedure  is  conducted  by  marching  along  the  frequency 
axis  in  fixed  steps,  filling  the  system  matrix,  and  calculating  its  determinant 
at  each  step.  A  change  in  the  sign  of  the  determinant  between  a  pair  of  adja¬ 
cent  frequencies  signifies  the  existence  of  a  root  between  these  frequencies.  A 

modified  secant  method  (ref  16)  is  then  invoked  to  accelerate  the  search  proce¬ 
dure.  The  root  is  determined  to  within  the  specified  tolerance  and  is  passed 
back  to  MOD-MAIN  along  with  current  values  in  the  system  matrix.  The  mode  shape 

coefficients  are  then  determined  by  a  call  to  the  coefficient  solver.  In  these 

routines  the  system  matrix  is  reduced  according  to  the  boundary  conditions 
imposed  on  the  left  end  and  the  coefficient  used  for  normalization.  The  new  set 
of  linear  equations  is  then  solved  for  the  unknown  mode  shape  coefficients  and 
control  is  passed  back  to  the  MOD-MAIN.  The  output  is  generated  by  calls  to  the 
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various  output  routines.  The  program  then  returns  to  the  top  of  the  loop,  the 
current  frequency  is  incremented,  and  the  process  is  repeated  for  the  next  root. 
The  looping  continues  until  all  of  the  requested  roots  and  coefficients  are 
found,  and  the  associated  output  is  generated. 

TEST  CASES 

This  method  was  evaluated  by  comparing  it  to  results  generated  by  an  FEM 
code.  The  software  chosen  was  ABAQUS  Release  4.7  developed  by  Hibbitt, 

Karlsson,  and  Sorenson  (ref  17).  The  beam  used  in  this  study  was  2540  millime¬ 
ters  (mm)  long  with  a  tapered  circular  cross  section  starting  at  a  left  end 
diameter  of  125  mm.  The  diametral  tapers  ranged  from  0.005  to  0.020. 

A  sketch  of  the  beam  and  its  two  modelling  representations  is  shown  in 
Figure  3.  The  finite  element  model  consists  of  53  nodes  and  52  two-noded 
linear  elements.  The  beam  elements  consider  the  effect  of  shear  but  not  rota¬ 
tory  inertia.  The  tapered  section  was  subdivided  into  13  prismatic  sections 
with  four  elements  in  each  section. 

The  USM  model  contains  three  segments  of  equal  lengths  whose  diameter  is 
the  average  diameter  of  that  section  of  the  beam.  Six  boundary  condition  com¬ 
binations  out  of  a  possible  32  were  evaluated  for  four  diametral  tapers.  The 
first  five  bending  frequencies  were  calculated  yielding  a  total  of  120  "data" 
points  for  each  model.  The  entire  120  runs  were  accomplished  in  one  afternoon. 

A  comparison  between  each  pair  of  values  indicates  that  at  most  a  6.7  per¬ 
cent  discrepancy  exists  between  the  results  of  both  models.  Most  calculated 
values  were  well  within  this  extreme.  In  general,  the  FEM  values  were  greater 
than  those  obtained  from  the  USM  for  the  first  three  modes.  However,  the 
reverse  is  true  for  modes  four  and  five.  These  comparisons  are  shown  graphi¬ 
cally  in  Figures  4  through  6  for  three  of  the  boundary  conditions.  In  these 
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graphs  the  frequency  is  plotted  against  diametral  taper  for  all  modes.  The 
upright  triangles  represent  the  values  determined  by  the  USM,  whereas  the 
inverted  triangles  are  results  calculated  by  the  FEM.  Overall  results  indicate 
that  the  effect  of  the  taper  is  greater  for  the  higher  modes.  For  example,  for 
the  fix-pin  beam,  the  fifth  frequency  is  reduced  from  970  to  820  Hz  as  the  taper 
is  increased  from  0.005  to  0.020,  a  15  percent  reduction.  The  primary  fre¬ 
quency,  however,  has  been  reduced  by  only  8  percent  from  60  to  55  Hz. 

In  Figures  7  through  9,  the  fourth  bending  mode  shape  of  each  model  is 
compared  for  the  extreme  taper  value  (0.020).  The  shapes  are  very  similar,  but 
the  magnitudes  of  the  USM  displacements  are  slightly  greater  than  those  from  the 
FEM.  The  greatest  discrepancy  of  11  percent  occurs  for  the  fix-guide  boundary 
condition  at  the  three-quarter  point  along  the  axis.  In  retrospect,  the  USM 
model  could  have  been  refined  by  either  changing  the  geometric  properties  of 
each  segment  or  by  adding  more  segments  to  better  approximate  the  taper. 

To  compare  the  differences  in  mode  shapes  as  a  function  of  taper,  consider 
the  results  plotted  in  Figure  10.  This  figure  contains  the  fourth  mode  shape 
generated  by  the  USM  for  the  fix-pin  boundary  condition  at  all  values  of  the 
taper.  As  is  expected,  the  beam  possessing  the  thinnest  cross  section  at  the 
right  boundary  has  the  greatest  deflection.  The  deflection  for  the  thinnest 
beam  is  about  15  percent  greater  than  for  the  heaviest  beam.  A  frequency  shift 
from  650  to  550  Hz  is  also  indicated. 

To  compare  the  differences  between  beam  modelling  POE,  consider  the  results 
shown  in  Figure  11.  In  this  figure,  the  mode  shapes  for  the  fourth  bend  mode 
and  the  steepest  taper  are  plotted  for  all  three  models  using  the  USM.  The 
boundary  conditions  are  fix-pin.  The  difference  between  the  shapes  is  hardly 
discernible  for  this  geometry.  The  frequency  shift  is  about  20  Hz,  the  highest 
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being  associated  with  the  Euler  model.  This  type  of  response  is  expected  since 
inclusion  of  additional  terms  in  the  POE  tends  to  retard  all  flexural  free 
vibration  frequencies. 

DISCUSSION  ANO  CONCLUSIONS 

This  work,  which  was  developed  in  conjunction  with  the  author's  disser¬ 
tation  (ref  15),  presents  an  accurate  and  efficient  method  to  determine  the  free 
vibration  frequencies  and  mode  shapes  of  nonprismatic  end-mounted  beams  on 
standard  boundaries.  The  term  Uniform  Segment  Method  has  been  coined  to 
distinguish  this  modelling  from  the  more  popular  finite  element  techniques.  A 
dedicated  computer  routine  has  been  developed  and  tested  against  results 
obtained  from  an  FEM  model  for  a  tapered  beam.  The  maximum  discrepancy  for  any 
frequency  calculation  was  7  percent,  whereas  the  maximum  difference  in  any  local 
value  of  the  mode  shape  was  11  percent.  Greater  accuracy  can  be  obtained  by  a 
modest  increase  in  the  number  of  USM  segments. 

Future  work  with  this  method  includes  linking  the  model  to  an  optimization 
package  so  that  discrepancies  between  USM-calculated  frequencies  and  a  set  of 
known  frequencies  (either  FEM-generated  or  experimentally  determined)  are  mini¬ 
mized.  The  optimization  variables  would  be  the  segmentation  spacing  and  the 
section  properties  of  the  segments.  Linearly  elastic  boundary  conditions  (in 
both  transverse  and  rotational  directions)  will  be  included  along  with  the 
standard  boundary  types.  Finally,  the  development  of  a  transient  analysis 
package  using  the  differentiable  USM  modes  will  permit  the  determination  of 
forced  responses  for  a  variety  of  displacement-dependent  loadings. 
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UNIFORM  SEGMENT  METHOD:  ANALYSIS  FLOWCHART 
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FIGURE  2.  UNIFORM  SEGMENT  METHOD:  ANALYSIS  FLOWCHART 


COMPARISON  OF  RESULTS:  FREQUENCY  RESPONSE 

UODES : 1 . 2 . 3  BEAM  BOUNDARY  COND I T I ONS :  F  I X  -  FREE 
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FIGURE  4A.  FREQUENCY  COMPARISON:  MODES  1-3;  FIX-FREE  BEAM 


COMPARISON  OF  RESULTS:  FREQUENCY  RESPONSE 
MODES:!, 2. 3  BEAM  BOUNDARY  COND I T I ONS : F I  X  -  GUIDE 


FIGURE  5A .  FREQUENCY  COMPARISON:  MODES  1-3;  FIX-GUIDE  BEAM 


COMPARISON  OF  RESULTS:  FREQUENCY  RESPONSE 

MODES : 4 , 5  BEAM  BOUNDARY  COND 1 T I ONS : FIX  -  GUIDE 
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FIGURE  58.  FREQUENCY  COMPARISON:  MOOES  4-5;  FIX-GUIDE  BEAM 


BEAM  BOUNDARY  CONDITIONS: 
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FIGURE  6A.  FREQUENCY  COMPARIS 


COMPARISON  OF  RESULTS:  FREQUENCY  RESPONSE 

MODES : 4 , 5  BEAM  BOUNDARY  CONDITIONS:  FIX  -  PH 
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FIGURE  6B .  FREQUENCY  COMPARISON:  MODES 


COMPARISON  OF  RESULTS:  MODE  SHAPE  RESPONSE 

TAPER:  .020  BOUNDARY  CONDITIONS:  FIX  -  FREE 
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FIGURE  7.  SHAPE  COMPARISON:  4 f h  MODE;  FIX-FREE  BEAM 


COMPARISON  OF  RESULTS:  MODE  SHAPE  RESPONSE 

TAPER:  .020  BOUNDARY  CONDITIONS:  FIX  -  GUIDE 
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FIGURE  8.  SHAPE  COMPARISON:  4 1 h  MODE;  FIX-GUIDE  BEAM 


COMPARISON  OF  RESULTS:  MODE  SHAPE  RESPONSE 
TAPER:  .020  BOUNDARY  CONDITIONS:  FIX  -  PIN 


25 


FIGURE  9.  SHAPE  COMPARISON:  4 1 h  MODE;  FIX-PIN  BEAM 


COMPARISON  OF  RESULTS:  MODE  SHAPE  RESPONSE 

MODEL:  USM  BOUNDARY  CONDITIONS:  FIX  -  PIN 


FIGURE  10.  SHAPE  COMPARISON  VARIOUS  TAPERS:  4th 


COMPARISON  OF  RESULTS:  MODE  SHAPE  RESPONSE 

MODEL:  USM  BOUNDARY  CONDITIONS:  FIX  -  PIN 

TAPER:  .020 
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